#!/usr/bin/python
"""Takes a bunch of sequences and returns desired number of sampled sequence.

Right now it contains two types of sampling; fermi dirac and boltzmann.
"""

# by Vineet Gupta <vineet@ural.wustl.edu>
# pdf.py 1-27-2007 13:52

# TO DO:
# 

import random

def fermi_dirac(sequences, percent):
    '''Returns list of sequences based on based on Fermi-Dirac distribution.
    sequences is a list of sequences out of which we need to generate
    samples. Percent is percent sequences out of given sequences as sample.'''
    
    samples_n = []
    samples = []
    
    #total_samples = len(sequences) * percent / 100
    total_samples = percent
    
    while total_samples > 0:
        samples_n = random.sample(sequences, 1)
        total_samples -= 1
        samples.append(samples_n[0])
        
    return samples

def boltzmann(freq_matrix, sequences, percent):
    '''Returns list of sequences based on based on Boltzmann distribution.
    matrix is matrix of real base frequency at various positions from which
    sample is to be generated. Actually it's just linear list here, but used as
    an array. sequences is a list of sequences out of which we need to generate
    samples. Percent is percent sequences out of given sequences as sample.'''
    
    bases = {'a':0, 'c':1, 'g':2, 't':3}
    samples = []
    samples_prob = {}
    total_prob = 0.0
    prob = 1.0
    
    #total_samples = len(sequences) * percent / 100
    total_samples = percent
        
    for seq in sequences:
        i = 0
        prob = 1.0
        for base in seq:
            prob *= freq_matrix[7*bases[base]+i]
            i += 1
        total_prob += prob
        samples_prob[seq] = prob
        
    while total_samples > 0:
        rand_num = random.random() * total_prob
        temp = 0.0
        
        for seq in samples_prob.keys():
            temp += samples_prob[seq]
            if temp > rand_num:
                samples.append(seq)
                total_samples -= 1
                break
        
    return samples
        
def main():
    pass
    
if __name__ == "__main__":
    main()
    